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We propose a phenomenological model for the polygonal hydraulic jumps discovered by Ellegaard 
et ai, based on the known flow structure for the type II hydraulic jumps with a "roller" (separation 
eddy) near the free surface in the jump region. The model consists of mass conservation and radial 
force balance between hydrostatic pressure and viscous stresses on the roller surface. In addition, we 
consider the azimuthal force balance, primarily between pressure and viscosity, but also including 
non-hydrostatic pressure contributions from surface tension in light of recent observations by Bush et 
al. The model can be analyzed by linearization around the circular state, resulting in a parameter 
relationship for nearly circular polygonal states. A truncated, but fully nonlinear version of the 
model can be solved analytically. This simpler model gives rise to polygonal shapes that are very 
similar to those observed in experiments, even though surface tension is neglected, and the condition 
for the existence of a polygon with N corners depends only on a single dimensionless number cj>. 
Finally, we include time-dependent terms in the model and study linear stability of the circular 
state. Instability occurs for sufficiently small Bond number and the most unstable wave length is 
expected to be roughly proportional to the width of the roller as in the Rayleigh-Plateau instability. 



PACS numbers: 47.35.-i,47.20.Ky,47.60.Kz, 68.03. Cd,ll. 30. Qc 



I. INTRODUCTION 

The term "hydraulic jump" refers to the sudden jump 
in fluid height, as for example observed in the outward 
spreading water layer in a kitchen sink, resulting from 
the impact of the water from the tap with the horizontal 
bottom of the sink [1-4]. Similar phenomena are seen 
in rivers with large tidal variation at the outlet and are 
known as "river bores" . River bores move up the rivers, 
whereas hydraulic jumps are stationary, either due to 
spatial inhomogeneities or to the geometric configuration 
of the flow. In the kitchen sink, the jump occurs on a 
more or less circular locus. Close to the point of impact 
the water level is very thin, but at a certain radius, rj(0), 
the level increases abruptly forming a circular jump. The 
water flow and thus the jump shape in a kitchen sink fluc- 
tuate, but by building a more symmetric experimental 
setup and/or using a more viscous fluid, one may obtain a 
completely stationary and axisymmetric flow, where the 
jump occurs on a surprisingly well-defined circle (FIG. 1). 
It has been shown [5-7] that circular hydraulic jumps can 
undergo a sequence of structural changes seen by varying 
the fluid height h a downstream of the jump. This can be 
achieved by inserting an adjustable circular weir at the 
rim of the circular impact plate. When increasing h a , 
the jump becomes steeper until, at a critical value of h Q , 
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the jump becomes unstable and loses its balance. Like 
a breaking wave, it creates a new stationary state with 
a wider jump region and a new flow structure in which 
the surface flow is reversed due to a separation vortex 
(referred to in the following as the "roller"). Following 
[7] , his new flow structure is called a type II jump to dis- 
tinguish it from the "ordinary" type I jump. The most 
remarkable observation about the type II state is that the 
jump typically loses its azimuthal symmetry and attains 
the shape of a regular, though not necessarily straight- 
edged, polygon (FIG. 2). 




FIG. 1. Circular hydraulic jump. A cylindrical jet of fluid 
impinges vertically on a horizontal glass plate and forms a 
circular hydraulic jump. The fluid used here is ethylene glycol 
(about 10 times the viscosity of water). 
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FIG. 2. Polygonal jump with eight corners, (a): Top view 
from an oblique angle, (b): View from below through a glass 
plate. The jump line becomes visible by diffracting light shone 
on the jump from above. The "belly" between two corners (or 
"necks") is where the roller is thickest. 

In this paper, we present a phenomenological model for 
the type II jump, which we believe to include the basic 
mechanisms of the flow rendering the circular state un- 
stable and giving the polygonal jump its peculiar shape. 
This allows us to obtain the qualitative dependence of 
the jump shape on the control parameters available, pri- 
marily the outer fluid height h Q , and the properties of the 
fluid, i.e., viscosity and surface tension. For experimental 
data, we refer to earlier work. Our own experiments are 
performed with the same setup as in [6] and is only used 
qualitatively as a visual guide (see also supplementary 
movies [20]). 

For the type I state, averaging theory with a variable 
profile gives at least qualitatively the right flow structure 
including the separation region on the bottom, outside 
the jump [8, 9]. For the type II state with separation 
at the surface, no such theory exists, and we therefore 
from the outset assume a flow structure derived from ex- 
periments. Our model extends earlier models proposed 
in [5, 6, 10], based on the competition between gravita- 
tional and viscous effects. In [5, 6] a radial force balance 
for the roller was used together with a "line tension" to 
derive polygonal shapes. In [10] it was shown that such a 
line tension is not necessary and that a more satisfactory 
model is obtained by taking into account the azimuthal 
flow. This is the approach used in the present paper. It 
was recently pointed out [11], however, that the instabil- 
ity forming the polygonal jumps seems to be driven, at 
least in part, by surface tension (see also Supplementary 
Materials [20]). Thus, surface tension is included in our 
model and indeed the instability forming the polygonal 
jumps seems to be closely related to the Rayleigh-Plateau 
instability for the "liquid cylinder" formed by the roller. 

In this paper, we treat both linear and nonlinear prop- 
erties of our model. For the linear properties, we can 
with some confidence include the effects of surface ten- 
sion, but for the fully developed polygonal states, this is 
very difficult because of the lack of precise data on the 
height profiles. Even without surface tension, our model 
does have polygonal states, and we compute their shapes 
assuming that they will be only slightly changed by sur- 



face tension effects. The layout of the paper is as follows: 
we first derive the model and point out the similarities 
and differences to earlier models. We then solve a non- 
linear model with zero surface tension analytically and 
show that it has many features that are expected from 
experiments. Finally, we compute the temporal stability 
of the circular state while including surface tension, and 
show that it has properties close to the Rayleigh-Plateau 
instability. 



II. DERIVATION OF THE MODEL 

In the experiment [5-7], the fluid falls from a nozzle 
mounted at some height, at constant volumetric flow rate, 
Q. We employ cylindrical coordinates, (r, 9, z), centered 
at the point, where the central axis of the cylindrical liq- 
uid jet impinges the horizontal plate. The height of the 
fluid surface above the plate is parametrized by h(r,0). 
As the fluid leaves the point of impact, h is small and the 
fluid spreads in a very thin layer (around 1 mm in our 
experiments) with supercritical flow speed (i.e. a speed 
larger than the speed of surface waves) . At a short dis- 
tance from the point of impact, the boundary layer in the 
film has become fully developed [2] and the fluid height is 
almost constant, i.e., h(r < rj, 9) f=a hi, as seen in FIG. 3. 
Due to the supercriticality, this inner flow is not affected 
by the transition from type I to type II jump, and the 
flow remains axially symmetric. The surface height of 
the type II jump abruptly increases at the jump radius 
Tj due to the presence of the roller on top of the thin fluid 
layer. The surface height increases for r,j < r until a cer- 
tain radius R where it settles at a roughly constant level, 
h(r > R,9) sa h . The height difference across the jump 
is Ah, and since h a ^> hi, we use the approximation 



Ah 



(1) 



throughout the paper. We further introduce the aspect 
ratio 



ho 
R 



(2) 



which is typically small, say of the order of 0.2. 

In FIG. 4 we show a visualization of the flow in a trian- 
gular jump created with ethylene glycol. It can be clearly 
seen that the roller fills the jump region from rj to R, 
where rj now has lost its azimuthal symmetry, whereas 
the outer edge R remains circular. We therefore define 
the normalized local roller width as one of the basic vari- 
ables of our model 



6(9) 



R- rj {9) 
R 



(0 < S < 1), 



(3) 



which varies with 9 to give the jump its characteristic 
polygon shapes. The dye dripped into the roller reveals 
an extremely slow exchange of the fluid in the roller 
(see also flow visualization in the Supplementary Ma- 
terial [20]). In fact the dye can be visible in the roller 
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FIG. 3. Measured height profiles of a pentagonal polygonal 
hydraulic jump in ethylene glycol. The solid curve represents 
the profile measured along the radial direction through a cor- 
ner. The dashed curve represents the profile along the radial 
direction midway between two corners ("belly"), where the 
edge of the jump is almost straight (see e.g. Fig. 2(b)). The 
height of the fluid layer was measured manually with a depth 
micrometer attached to a translation table. Reproduced with 
permission from [16]. The flow rate is Q = 40 ml s _1 . 
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FIG. 4. (Color online) Flow visualization of the roller vortex 
in the case of a triangular jump in ethylene glycol. The jump 
is seen from below through a glass plate, where the impinging 
jet is visible as the black center region and the jump line as 
the black line surrounding it. Red dye is injected into the 
roller flow by letting droplets of the dye fall from above into 
the vortex. It is seen that the roller structure extends from 
the jump line to the outer radius. In the corner region, fluid is 
expelled in clearly visible radial jets. Note also the fine white 
line between the roller vortices, which indicate that the vor- 
tices are actually disconnected structures. For a visualization 
of the flow see also Supplementary Material [20]. 



for several minutes if it is left undisturbed. From this 
we conclude that the main part of the fluid in the thin 
layer going under the roller leaves the jump region with- 
out ever entering the roller. Downstream of the jump 
region, the flow speed is subcritical (less than the speed 
of surface waves), and the flow is again purely radial, but 
it now has an azimuthal dependence: at the corners of 
the polygon, strong radial jets are observed. Within the 
jump region, an azimuthal transport must therefore ex- 
ist in the jump region from the sides to the corners of 
the polygon. Indeed a small, spiraling azimuthal flow is 
observed in the roller from its "belly" out to the "neck" . 
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FIG. 5. Measured phase diagram, reproduced with permis- 
sion from [6]. Polygonal hydraulic jumps are observed using 
ethylene-glycol, which has a viscosity about 10 times that 
of water. Polygons appear in the parameter regime between 
the circular type-I state (marked as T) and the closed state 
(i.e. no jump, marked as 'C') when the jet flux Q and the 
outer height h are varied. The height of the nozzle about 
the surface is set to 2 cm. The number of corners is found 
to be less sensitive to the nozzle height [6]. Polygonal states 
with the number of corners N = 1, 2, . . . , 8 are found in this 
experiment. 



From FIG. 4 the jets seem to consist primarily of liquid 
going right through the jump, but they also carry with 
them the dyed fluid transported to the corners inside the 
roller. The jets at the corners of an iV-gon are actually 
so strong that they apparently break up the roller into 
N distinct rollers, touching and interacting in a rather 
complicated manner in the corners [6] . This effect is not 
included in our model, but since it occurs only in nar- 
row regions near the corners, we believe that our model 
can still give useful results. A phase diagram has been 
presented in [5] in terms of three parameters: the flow 
rate Q, the outer fluid height h a and the height of the 
nozzle. The dependence on the first two parameters is 
reproduced in FIG. 5. When increasing h a , the tran- 
sition to the type II structure leads to a polygon with 
many corners (up to N = 13 have been observed). In- 
creasing h further, the number of corners decreases and 
so does the average jump radius, (2n)^ 1 <f r,j{9)d9 1 un- 
til at last the jump closes on the vertical jet from the 
nozzle. Throughout this regime the "outer" radius R 
is surprisingly constant, while the inner border defining 
the polygon changes. In fact, this polygon shape shows 
considerable hysteresis, and a whole series of bordering 
polygon shapes can be stable at the same flow rates. An 
extra corner can be created by pulling out a polygon side 
with a needle, and similarly, a corner can be removed 
by pushing a corner inward. For sufficiently small dis- 
turbances the jump simply deforms - just as if a rubber 
band was pulled - and bounces back to its original shape 
after being released. 
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FIG. 6. (Color online) A three-dimensional view of the in- 
finitesimal wedge [8, 8 + d8] for the hydraulic jump. The gray 
shaded region is a part of the roller eddy near the surface just 
after the jump, and represents the control volume we consider 
for the mass and force balances. Areas of the side, bottom, 
rear faces, as well as the free-surface portion of the volume 
are labeled as A t (9), dAt,(6), dA r , and dAf(8), respectively. 
The bold blue arrows indicate the flow direction. 



A. Momentum Balance for the Roller 

In our phenomenological model, we focus on the roller 
which we consider as a separate object interacting with 
the fluid beneath and behind it. We look at a control 
volume defined by the section of the roller inside an in- 
finitesimal wedge [8, 9 + dd], as illustrated in FIG. 6, and 
label the areas of the bottom and the rear faces of the 
control volume as dAb(8) and dA r , respectively. The lat- 
eral areas are denoted as A t (8). Let further be the 
stress tensor, let V and S be the volume and the bound- 
ary, respectively, of the roller slice, and let n be the unit 
normal of S. Then the momentum equations are 

J ^ipa)dV + J pa{u-h)dS 

= - f p(ndS)+ f ry(n,-dS)+ / pfdV. (4) 
Js Js Jv 

Initially, we shall consider stationary states and thus the 
first term on the left hand side is zero. Also, the last 
term on the right vanishes, since the only body force is 
gravity, which is taken into account by assuming that the 
pressure p is hydrostatic. 

The inner flow in FIG. 6 is purely radial and axially 
symmetric. When it reaches the roller, most of it passes 
underneath it. When it has passed the roller the flow is 
still radial, but has acquired an angular dependence. In 
between, in the region of the roller, the flow thus acquires 
a small tangential component, which we assume to run 
inside the roller. Our first task is then to write down the 
continuity equation for the flow. Then we shall evaluate 
the "force" terms in (4) for the radial and tangential di- 
rections, respectively. In the following, we use cylindrical 
coordinates and denote the respective vector components 
of dependent variables by indices r, 8 and z. 



B. Mass Conservation 

FIG. 7 (a) shows the top view of the entire roller, and 
Fig. 7 (b) displays an angular section dd of the roller, 
together with the mass flux entering and leaving the con- 
trol volume (corresponding to the gray shaded volume in 
FIG. 6). Since the flow before the jump is purely radial 
and independent of 8, as evidenced in FIG. 3, the radial 
flow into the volume must be ^d8. The fluid going out of 
the volume in the radial direction is denoted by q r {9)dd. 
The fluid going into the roller at the cross sectional wall 
at 8 is qe(8), and the corresponding flow out of the roller 
at 8 + dd is qg (8 + d8) . Mass conservation in the region 
then demands 

~dd + qg(8) = q r {8)d8 + q e (8 + dd). 
Ztt 

We define the normalized radial and azimuthal flux per 
azimuthal angle £ r = q r /Q and £g = qg/Q, respectively. 
The continuity equation is now written as 

By integrating from 8 = to 2tt and using the periodic 
condition for we obtain a constraint for the total flux 
out of the jump region. 

I £ r (6)d9 = 1. (6) 




FIG. 7. (Color online) (a) Schematic top view of a pentag- 
onal jump. The gray shaded wedge represents the control 
volume with infinitesimal angle dd, also seen in Fig. 6. (b) 
Enlargement of the wedge-shaped control volume. Arrows 
(blue) represent the radial and tangential mass flux q r and 
qe, i.e. the flux in and out of the control volume. 



C. Radial Force Balance 

We now inspect FIG. 8, the side view of FIG. 6, and 
construct a force balance equation in the radial direction 
for the roller which we view as a stationary body of liquid. 
It experiences a force directed radially inward due to the 
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difference in fluid height at its inner and outer rim. The 
shear force also acts on the bottom of the roller, caused 
by the velocity gradient of the radially outward flow along 
its bottom boundary which we estimate at z — hi. We 
assume that these forces balance, and that surface ten- 
sion plays a secondary role in the radial direction. Under 
this assumption, we have 



dF? + dFlf = 0, 



(7) 



where dF^ and dF^ are, respectively, the hydrostatic 
pressure force and the radial component of the viscous 
friction force per unit azimuthal angle on the roller slice 
in FIG. 6. 
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FIG. 8. (Color online) Side view of the wedge in Fig. 6. The 
height inside and outside of the jump are hi and h , respec- 
tively. We treat the roller eddy (shaded region) as a separate 
body of fluid from the main stream, and assume the force 
balance of the hydrostatic force dF^ and the shear force dF^ 
in the radial direction (forces indicated by bold red arrows). 



1. Hydrostatic pressure force. 

The pressure force, dF h , is obtained by integrating the 
pressure over the area of its outer rim assumed to be at 
r = R and extending from z — hi to h a (see Fig. 8). 
Further, assuming the pressure to be hydrostatic, and 
omitting the ambient pressure which cancels, the force 
becomes 



dFl L = - 



hi 



pg{h — z) dz Rd6 



-pgR^de, 



(8) 



using (1) where p is the fluid density and g is the gravi- 
tational acceleration. 



2. Shear. 

If the fluid is Newtonian with the dynamic viscosity /z, 
the outward force on the roller dF^ shown in Fig. 8 is 



obtained by integrating the shear r r 
the bottom area of the roller: 



du r 



r drd9, 



over 



(9) 



where u r is the radial velocity component. To estimate 
this integral, we assume that the thin radial flow contin- 
ues forward under the roller in a height of the order h* , 
an appropriate characteristic height near the roller. The 
velocity in this region is related to the radial flux q r by 
q r = u*r*h* where r* and u* denote some characteristic 
radius and radial velocity, respectively. We also assume 
that the shear can be expressed as 

du r 
dz 



-ci 



ciq r 



z—hi 



h* 



with a numerical factor c\ > of order unity. Then, we 
have 



(R 2 -r 2 )d9. (10) 



For simplicity, we choose r* = |(r^ + R) and h* = h t . 
Then the shear force becomes 

dF? « °^{R - rj )de = c ^ R ^ r Sde. (11) 



3. Radial balance. 

Substituting the estimated forces (8) and (11) into 
Eq. (7), we obtain the radial force balance equation 

with the nondimensional parameter 

ghfhl 



(12) 



III 



(13) 



2c\vQ 

where v is the kinematic viscosity \ij p. This equation 
was first derived in [5]. It shows that the radial flux 
increases at the corners where S decreases. For a very 
thin corner, the flux would increase to infinity, and al- 
though our model remains meaningful in this limit, we 
do not expect it to represent the physics correctly, since 
the detailed structure of the corner flow could play an 
important role. Lacking such a more detailed model, we 
shall thus keep in mind that polygons with very thin or 
vanishing rollers in the corners are probably somewhat 
idealized. 

In the circular state the normalized radial flux is £ r = 
1/2-7T so the dimensionless roller has width So where 

S = 2ttIIi. (14) 

Combining (12) with the mass conservation (5), we 
obtain 



— = 5i 

2tt S 

whereas prime denotes the derivative with respect to 9. 
The constraint (6), along with Eq. (12), becomes 

f d9 _ 1 



(15) 



(16) 



G 



D. Azimuthal Force Balance 

We now discuss the azimuthal transport in the con- 
trol volume shown in FIG. 7. Since the flows in the az- 
imuthal direction are slow, we neglect the kinetic term 
arising from the second term on the left hand side of (4) , 
and consider hydrostatic pressure, viscous, and surface 
tension terms. 



Considering the directions of the forces acting on the 
sides, the force difference between the two sides is 

dF ° = ~Te Feh(d)dd * -^1T 6/{d)d9 ' (22) 

where prime again denotes the derivative with respect to 
9. Due to the minus sign, the force acts in the direction 
from a belly to a corner of the roller. 



1. Kinetic term. 

The second term on the left hand side of (4) is the 
momentum flux through the boundary of the roller slice. 
We consider the flux for the vertical sides of the cross sec- 
tions at 9 and 6 + d9 in FIG. 6. Using the approximation 
h 3> hi, the areas of these sides are given by 

A t (9) a RAh8(6) a Rh 6(9). (17) 

We consider the average azimuthal flow velocity over such 
a surface ug — qg /A t . Then, the kinetic term is estimated 
as 

dF e km = [ P u^ 2 A t ] e e +de a puo- 2 Rh 5(9)d9 (18) 

Since the mean velocities in the roller always remain 
small, except, perhaps inside a very thin corner (where 
our model breaks down anyway, since the roller breaks) 
we shall neglect this term in the following. 



3. Shear. 

The contributions to the azimuthal component of the 
second term on the right hand side of Eq. (4) are the 
shear forces between the azimuthal flow inside the roller 
and the radial flow underneath and behind it. Just as 
for the radial force balance we must make some estimate 
of the velocity gradients involved. The outer rim of the 
roller has area 

dA r « RAhdO « Rh d9 (23) 

on which we estimate dug / dr\ r=R w —C2Ue/(R—rj) since 
uq ~ for r > R. Here, c 2 is a positive numerical factors 
of order unity, and ug is the average azimuthal flow ve- 
locity over A t . The areas of the sides are given by (17). 
Thus, 

W=qe/A t ^qe/(Rh 5). (24) 



2. Hydrostatic pressure. 



Likewise, the area on the bottom of the roller is 



For the azimuthal component of the first term on the 
right hand side of Eq. (4), we assume the pressure in the 
roller to be hydrostatic. Then, a net azimuthal pressure 
force on the roller slice arises due to the difference in area 
of the two vertical sides at 9 and 9 + d9 shown in Fig. 6. 
The force on one side is estimated as in Eq. (8). For this, 
we assume a simple surface profile which is linear within 
the cross-section 9 = const., i.e., 

Hr ' 6) Khl + Wf) {r ~ R + Rm) (19) 

where h = hi at r = rj(9) = R(l — 6(6)), and h = h at 
r = R. With h 3> h l we get 



h(r,9)^—{r-R(l-5(9))} 
with the aspect ratio defined in (2). We then get 



(20) 



nR fh(r,8) 

?h (a \ _ / dr 



FS(9)= / dr / pg(h(r,6)-z)dz 

Jr{x-s{6)) Jo 

pgRh 2 







5(6). 



(21) 



dA h 



r])d9 



(25) 



on which dug/dz\ z=h . « c 3 ug/h* sa c 3 ug/hi since ug = 
on z = 0. Here, we have used the height hi rather than 
Ah to estimate the shear between the roller and the main 
flow, and C3 is another positive numerical factor of order 
unity. 

Note that the radial velocity gradient at the rim is 
measured from the interior to the exterior of the roller, 
whereas the velocity gradient on the bottom is in the 
opposite direction. Thus, we reverse the sign of the latter 
in order to add both contributions for an estimate of the 
shear force on the control volume: 



-pug 

' Rh S 
pQ£,e 



c 2 , c 3 R - r } 



R-r 



■Rh 



d6 



cih c 3 RS(R + rj) 
~5~ + 2h~i 
c 2 c 3 R(2-8) 



d6 



RS 2 1hih 
Here, the relation (24) is again used. 
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4- Surface tension. 

Surface tension acts by changing the pressure inside 
the roller depending on the local curvature. This does not 
to a first approximation change the radial force balance, 
so we shall concentrate on the effects in the azimuthal 
direction. A somewhat similar approach in a different 
context can be found in [12]. To model this, we shall 
include the Laplace pressure difference across the surface 
z = h(r,8): 



Ap = a 



1 

R~i 



1 



(26) 



where R\ and R 2 are the two principal curvatures; we 
take their signs to be positive if the center of curvature 
is on the "inside" direction from the surface (i.e. located 
under the surface). Taking the atmospheric pressure to 
be zero everywhere, the relation (26) describes the fluid 
pressure just beneath the surface. Thus, the total pres- 
sure inside the fluid becomes 



p(r, 9, z) = pg(h(r, 9) - z) + Ap. 



(27) 



To use these expressions we need information about 
h(r,9), in the circular state and in the polygonal state 
— and in the family of states between them. Of course, 
we do not have this information, so, in keeping with our 
simple model, which expresses the geometry of the jump 
entirely in terms of the local width of the roller, we shall 
replace the two radii of curvature with averages over the 
cross-section of the surface. 

In what follows, we consider small perturbations, i.e. 
nearly circular jumps. One principal radius of curvature 
i?i quantifies the curvature roughly in the (r, z)-plane, 
i.e., the cross section of the jump. This is the "danger- 
ous" one, which may lead to instability. The other one, 
i?2 , is defined in a plane orthogonal to R\ , also including 
the surface normal. Since surface inclinations are small, 
i?2 can be approximated by the radius of curvature of the 
jump shape seen from above (cf. Fig. 7). For slow varia- 
tions of the plane curve R5(9) the radius of curvature we 
get 



1 

~R~2 



S"(9), 



(28) 



where primes denote differentiation with respect to 9. 

To estimate R\ we have to know how the shape distorts 
when a corner is created. Since it is observed that addi- 
tion of surfactant destroys the polygons (see also Supple- 
mental Material [20]), it seems likely that surface tension 
can make the circular states unstable through a Rayleigh- 
Plateau-like instability (however, it cannot explain the 
observation of circular type-II states that are sometimes 
stable.). This means that the curvature should increase 
when S is decreased as it does near a corner. This is 
indeed observed in experiments as seen in Fig. 3, where 
the height profile (on a radial path) of a polygon-state 



has been measured with a conducting needle [16], both 
through an edge and a corner. It is clear that the cur- 
vature is larger (i.e. the radius of curvature is smaller) 
through the corner. Thus, we assume 



R[(S) > 



(29) 



where prime denotes differentiation with respect to 5. 

When the height Ah and the width R5 are similar, we 
expect that R\ would have a value close to these, like for 
a cylindrical surface, but in contrast to the classical cylin- 
der case of Plateau and Rayleigh, the curvature cannot 
keep increasing when the roller shrinks since the external 
height is fixed. Thus, we expect the curvature to saturate 
when the width of the roller becomes much smaller than 
Ah. For future use, we shall use the symbol p\ for the 
normalized cross-sectional curvature: p\ = R\/R. The 
Laplace pressure term in (27) now becomes 



Ap = 



1 



R \ Pl {5) 



-5" 



(30) 



To get the force, we now have to integrate Ap ■ n over 
the bounding surfaces, where n is the local outward nor- 
mal, and project the result along 9. This will give the 
force from the roller on the surroundings. To get the 
force on the roller we have to put in a minus-sign. The 
rear surface does not contribute anything since it is par- 
allel to 9. Treating Ap as constant in the cross-section 
9 = const, yields 



F v . 



-A t A P Q 



(31) 



where the areas of the sides A t are given by (17). 

To estimate the free surface contribution we need the 
area of the free surface and the ^-component of the out- 
ward unit normal vector. As above, we use the linear 
height profile (20), where 



and 



dh _ Ah 



_ dh _ 5'{Q)Ah f r 
,e= d9~ 6(9) 2 V~R 



(32) 



(33) 



where, for brevity, we use a compact notation for partial 
derivatives throughout this section. The area of the free 
surface is 



dAf k Rd9^(RS) 2 + (Ah) 2 , 
and the unit normal vector is 



n = V j -h,e /r 



(34) 



(35) 



with a suitable normalization factor N. Assuming that 
S'(9) <C 1 so that \h,e /r\ -C h, r and \h,g jr\ -C 1, we have 



N : 



1 



R8 



v^Mp + (Rsy 



(36) 
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The ^-component of the unit normal becomes 

S h,g 



Nh, e 
R 



(RS) 2 ) 

and we obtain (based on the assumption Ah <C RS) 



(37) 



dAf hg « — RS h 
Using the average value for h,$, 
1 



dO. 



RS(6) 



R(1-S(S)) 



we obtain 



dAf hg 



1 



2 6(0) y '' 



R h a 5' did. 



(38) 



(39) 



(40) 



The force contribution from the free surface is this prod- 
uct times Ap. Thus, 



dF tt , 



—dAf fig Ap 
uh , ( 1 



Pi(S) 



s" de. 



(41) 



The total surface tension force in the azimuthal direction 
can now be computed from (31) and (41). Note that we 
need the forces on the roller, so the direction of the force 
is into the roller, i.e., in the direction 9 at 9 and in the 
direction —9 at 6 + d9. Thus, 



dFS 



-F' 



-h D crS 



t ^(9)d9 + dF !: 
S'p[(S) 



pi 



+ S"'(0) 



(42) 



It is seen that if p[(So) > (like in the case of a cylinder 
squeezed uniformly), the first term gives a contribution 
similar to the hydrostatic one (22), but with opposite 
sign. This seems to be the case for the polygonal hy- 
draulic jump shown in FIG. 3, where the curvature in 
the corner is larger than at the edge. 

For later use, we shall introduce the shape parameter 



B 



Sop'xjSo) 
(Pi(6 )) 2 



(43) 



For a cylindrical roller with p(S) « S we would then get 
B « Sq 1 . The ratio of the prefactor of the gravitational 
term (22) to the surface tension term (without S) is the 
dimensionless number 



pgRh 2 /2 Rh 



Bo 2 
3a 



(44) 



where L 



3a h /2 3l 2 
[af(pg)] 1 ^ 2 is the capillary length and where 



Bo 



ho 



(45) 



is the Bond number based on h n 



5. Azimuthal force balance. 

Using these estimates, we write down the azimuthal 
component of Eq. (4): 



dFS 



dF? + dFS 



0, 



(46) 



which can be written in dimensionless form as 

5' 



S'{9) = 



II2 

S 2 



n, 1 



ft 



+ 3aBo~ 2 c5 



Pi 



S'{9) + 5'"{9)\, (47) 



where we have defined the dimensionless numbers 
2c 2 vQ , _ 2c 3 vQ 



IF 



and II 



ghihl 



(48) 



The two coupled Equations (47) and (15) close our 
model for 5 and ft . An overview of the numerous symbols 
used in our model is provided in Table I in Sec. B. 



E. Linear analysis 

We investigate the existence of nearly circular poly- 
gons, i.e. for small oscillations around the circular solu- 
tion 5 = So = 27TIT! and ft = 0. Inserting 5 = So + e#i(0) 
and ft = eft(#) into (47) and (15), and by using B in 
(43), we find the first order perturbation equations to be 



(1 - 3aBo~ 2 B) 6[ = -ft 



+ 3aBo" z 5o S'" 
zirdf) 



!+n s( i-i j0 ) 

-2; 



(49) 
(50) 



Differentiating (49) with respect to 9 and substituting ft 
from (50), we obtain 



(1 - 3aBo" 2 B) S'l = 



1 



2nS Q 



II9 1 

^ + n 3 (i--5 ) 



S 3aBo~ 2 S'l' 



Si 

(51) 



For perturbations with Si = sin(fc(# — #o))i we obtain the 
characteristic equation: 



S k 4 + k 2 



Bo 2 
3 a 



— B 



Bo 2 
6irado 



This equation has the solution 



k 2 = 



25 



6a 



= 0. (52) 



(53) 
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where we define 



G 



D 



Bo 2 
3 a 



2Bo 2 
3a7r 



II2 

S 2 o 



n, 1 



~-6 



(54) 



When the wave number k is an integer, the solution is 
periodic on [0, 2Tr/k] and corresponds to a polygon with 
k corners. In the present case, k 2 needs to be positive 
for the solution to make sense. In the limit of negligible 
surface tension where a — > (i.e. Bo — > 00), we simply 
get 



k = ±A 



1 



II2 

L*2 



n. 



1 - -S 
2 



(55) 



III. SOLVABLE NONLINEAR MODEL 

To solve the static nonlinear model, Eqs. (47)-(15), we 
would need to know pi(S). At present we do not know 
the shapes well enough to estimate this, and we shall 
restrict our attention to the case where surface tension is 
absent; despite this simplification, we shall see that we 
still obtain meaningful results which allow us to explain 
the structure of the phase diagram. 

It is instructive to start with a further simplification. 
By letting il 3 — > 0, we neglect the term due to the shear 
against the bottom of the roller. By ignoring this regular 
perturbation term, the model can be solved analytically, 
which enables further analysis, as discussed in Sec. HID, 
where we investigate the effect of the 1I3 term. We obtain 



0V 



y' = i- 



1 



where we have introduced the rescaled variables 

5 
So 

y{0) = 27T&, 



x{6) = 



and the non-dimensional parameter 



(27r) 2 n? /2 n 2 - 1/2 



n 2 hlhjg 2 R 

V 2 Q 2 yJc\c 2 



(56) 
(57) 

(58) 
(59) 

(60) 



We now solve (56) for y, differentiate it once with respect 
to 9, and eliminate y' in (57), i.e., to obtain the single 
(second-order) equation 



1 



1. 



(61) 



By substituting X — x 3 /3, this may be formulated as an 



equation of motion for a 
force field, 

(b 2 X" 



mass 



dV_ 
~d~X'' 



in the conservative 



(62) 



where the potential function is given by 

V = X--(ZX) 2 ' 5 = — - 
2 V ' 3 



We introduce the rescaled angle G 
once in to obtain 



(63) 



and integrate 



(64) 



where C is a constant of integration corresponding to the 
"energy" of the mass c/> 2 . Solving the equation for dX/dO 
and transforming back to x, we obtain 



dx 
dO 



± 



V2(C - V(x)) 



This equation has the implicit solution 

e(x b ) = e(x a ) ± f 

J X 



y/2(C-V(x))' 



(65) 



(66) 



With (58) we recover the expression in terms of the orig- 
inal scaling 



S/S 



x 2 dx 



8(6)=9(S a )±cb (67) 
JSa/So v / 2(C~V(x)) 

and similarly, by substitution of (65) into (56), the cor- 
responding azimuthal flux is given by 

ie = T^2{C-V{8ISo)). (68) 



A. Phase Diagram 

We now consider the parameter space for existence of 
polygons. Such solutions may be considered "bound" 
states of the potential V. As shown in FIG. 9, such a 
state exists if and only if — i < C < 0. When C is within 
the range, there are two roots x min and a; max for the equa- 
tion V(x) = C, satisfying < x min < x max < |. They 
serve as turning points for a trajectory, and describe the 
minimum (corner) and the maximum (belly) of the jump 
width, respectively. Without loss of generality we may 
choose x = x min when 9 = = 0, and then a trajectory 
can be computed by integrating (65) with the plus sign 
from 6 = to 6 = T/2, where T(C) denotes the period 
of oscillation in terms of Q, and extending that part of 
the solution using symmetry. 

For a polygon with N corners a trajectory must oscil- 
late in the potential N times, and this must result in an 
increase by 2tt in terms of 9 = <pQ. Thus, a solution for 
iV-polygon must satisfy a commensurability condition 



<j)N 



2tt 

nc~y 



(69) 
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x,X 



FIG. 9. The potential function V as a function of x — 5 /So 
and X — x 3 /3, respectively. A periodic solution exists only 
when — | < C < 0. The solution corresponds to a polygon 
with N corners if the period is 2n/N in terms of 6. 



with the normalized half period 
T(C) 



«(C) 



x dx 



l(C ) y/2(C-V(x)) 



(70) 



When C — — i the equilibrium solution x mia = x max = 
1, corresponding to the circular jump 6 = 8q, is obtained. 
When C is only slightly larger than -| wc approximate 
the potential near the local minimum as V(x) ~ — g + 
\{x — l) 2 . If we express x min = 1 — e, then ir max = 1 + e 
and the "energy" is approximated by C = V(a; min ) = 



V(x max ) ~ - 
and we find 



2 



2 fc ■ 



lim 



Then, 2(C - V(a;)) ~ e 2 - (a - 1) 



1+6 



(71) 



for which the condition (69) becomes <j> = 1/N. Nat- 
urally this agrees with the condition (55) found in the 
linear analysis when II3 = 0. 

From C = — g to C = 0, the oscillation amplitude 
grows until, at C — 0, x min = and x max = |. Corre- 
spondingly, the integral T(C)/2 decreases monotonically 
in C until it reaches 



T(0) 



3/2 



x dx 



yJ~2V(x) 



= 3. 



(72) 



In this limit the condition (69) becomes <f> = ir/(3N). 

To summarize, the condition (69) for the existence of 
a polygon solution becomes 



K 

N' 



(73) 



where K = 2tt/T e [1,tt/3] « [1,1.0472]. if = 1 corre- 
sponds to nearly circular and if = 7r/3 to spiky polygons. 

We also need to ensure < 6 < 1. Since x„, in > 0, 
S > is guaranteed. Thus, the remaining condition is 
5 < 1 for all 9 which is equivalent to 



N(f) 
it/3 







spiky solutions (S) 


\ 

\% 


circular solutions (C) 






0.4 


0.8 



5 < l/a; 



(74) 



5 



FIG. 10. Polygonal jumps with N corners exist in the pa- 
rameter range shown above (truncated nonlinear model), in 
the truncated nonlinear model. The jumps bifurcate from 
the line cj> = 1/N for the circular jumps (C). When the pa- 
rameters approach the line (f> — n/(3N), the roller thickness 
approaches S — in the corners at some 9, indicating a rup- 
ture; solutions appear spiky near the line (S). On the border 
So = l/i m «, the roller thickness reaches S = 1 so the solution 
becomes unphysical (UP). 



On this border the roller thickness is predicted to equal 
the jump radius. This unphysical behavior is likely to 
be due to truncation of the terms in this approximation. 
The border curve is computed in FIG. 10 together with 
the lines <j> = l/N and ir/(3N). 

Since £ max is determined purely by the choice of C, 
the border (74) is identical for all N but only scaled and 
shifted in the (^-direction. Apart from this border, an A^- 
polygon exists in the horizontal strip n/(3N) > <f> > 1/N. 
It is natural to ask for which A^ we find an overlap of 
these bands corresponding to multistability of polygons 
as observed experimentally. Consider two neighboring 
strips with A^ and N + 1. The condition for overlap of 
the strips is </>^ in < 0^+ l5 or 1/iV < n/{3(N + 1)}. This 
leads to a condition N > 3/(tt - 3) « 21.2, oriV> 22. 
While the model thus in principle predicts overlapping 
regions and resulting hystereses of polygonal jumps with 
different N, it is an experimental fact that overlap exists 
for much lower N. 



B. Solutions and Shapes 

Some typical solutions are displayed in FIG. 11 in 
terms of the variables 6(9) and q r (9) for the symmetries 
N = 1,2, 3, 4, which were obtained by solving Eqs. (56)- 
(57) numerically. In Fig. 11, each row shows the transi- 
tion from nearly circular to spiky solutions - this may be 
achieved by tuning N<j> between 1 and tt/3. Notice that 
a solution only exists for a special value of the energy 
C, and is determined via Eq. (69) by the choice of the 
symmetry A^ and "mass" <fr 2 . The strict monotonicity of 
T = T(C) in C (see above) means that this mapping is 
one-to-one. It is meaningful to start integration in a cor- 
ner where x' = at x min or x max , respectively. The initial 



11 



condition x mi „ is then given by this C, (see Eq. (65) which 
yields V^(a; mInjlnajc ) = C for x' = 0). In other words, to ob- 
tain the desired symmetry N with mass <fi 2 , the initial 
condition x min must be chosen accordingly. 




FIG. 11. (Color online) Typical polygon solutions are shown 
for the truncated model in Eqs. (56)- (57). The polygon shape 
(jump line) is given by 5(9) (solid red) and the radial flux by 
£r(0) (blue dashed), which is computed via the radial balance 
Eq. (12). Numbers refer to the symmetry N and increase 
down the columns. From left to right, the transition from 
nearly circular to spiky is shown for each symmetry (V, 'm' 
and 'h' refer to center, middle or homoclinic orbit, respec- 
tively), corresponding to the transition discussed in Fig. 10. 



C. Comparison with experimental observations 

We now represent the shaded region in FIG. 10 in 
terms of the physical entities. Our aim is to compare the 
predicted region with the measurement on the (h , Q)- 
plane in FIG. 5. The two parameters in the model can 
be written as 



S = irghfhl/(civQ) 

i 2 2di3l4;/ 3/2 1/2 2/^,2 



(75) 
(76) 



We notice that both 5q and <p appear to depend on h 
and Q only through the combination h 2 Q /Q. If so, then we 
are not able to map the bifurcation curves from FIG. 10 
one-to-one to the phase diagram curves in FIG. 5. 

However, the jump radius R actually depends on 
other parameters as described before. For an esti- 
mate of this dependence, we use the relation R = 
C4 : Q 5 / 8 is~ 3 / 8 g~ 1 / 8 (c 4 = const.) proposed in [15], since 
it is based on experimental data and leads to a sim- 
ple estimate for the mapping. This converts Eq. (76) 
to ~ 7r 2 C4 5 15 / 8 ^^/(c? /2 C 2 /2 ^ 19 / 8 g 11 / 8 ), and the Q- 
dependence of <j> is corrected. Solving this with (75), we 
obtain 



cic\v 8 (j) 8 

7T 5 c5(7 4 /l?5 U 



Q 



4gp 3 h 8 (j) s 

c l c 4°0 



1/10 



(77) 



1/5 



Thus we are able to map (<5o,i 
parameters h a and Q. 



(78) 

back to the physical 
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The spiky polygons go along with a singular behavior 
of the radial outward flux q r (0) (blue). This behavior 
is caused by 5 —> in Eq. (12), and it is likely that this 
behavior disappears by inclusion of some of the neglected 
terms. However, we note that although this behavior is 
not physical, it is seen that the polygons exhibit very 
strong jets in the corners in the outward radial direction, 
as is seen in the flow visualization in Fig. 4. Notice also 
that the wave number has the right dependence on h , 
i.e. larger h Q leads to fewer corners. Thus the increase 
in the number of corners with increasing Q — as seen in 
FIG. 5 — is reproduced. 



FIG. 12. Bifurcation diagram in physical parameters (h , Q). 
The regions in which polygonal solutions with each iV up 
to TV = 8 are shown. The h -Q relationship for the bor- 
der curves, especially in the enlarged inset, is qualitatively 
similar to the measured FIG. 5 despite quantitative disagree- 
ment. The bands are predicted to be thinner and do not 
overlap in the model, unless JV > 22, unlike in experimental 
observations. 

When only Sq is varied while <f> and other parameters 
on the right hand sides of (77,78) are fixed, these relations 
predict h oc <5 n / 10 and Q ex. 5 16 / 5 : hence Q oc h^^ 11 . 
Thus, the top and bottom border lines in FIG. 10 are 
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mapped to these curves. Points on the left border line 
Sq = diverge in the physical plane, and the remaining 
border Sq = l/x max is mapped to a curve close to the 
origin of the physical plane. For each N = 1, 2, . . . , 8 we 
have mapped the region in FIG. 10 in this way, and the 
resulting physical parameter plane is shown in FIG. 12. 
Here, the typical set of parameter values in Appendix A 
are used for physical parameters other than h Q and Q. 



D. Nonlinear model with II3 > 

An analysis of Eqs. (15) and (47) where H3 > is also 
possible, by numerical means, and partially also by using 
more intricate analytic arguments based on a Hamilto- 
nian formalism which we only mention here for the sake 
of brevity. 

To numerically solve the static equations including 
II3 > 0, a scheme similar to the one described above must 
be adopted. For this case, we have no analytic means 
of determining the period T(C) in the same manner as 
above. However, one may integrate the equation system 
with an event solver: if we start integration from the 
corner, a half period is given when S' = (corresponding 
to the "belly" of the polygon), which defines the condi- 
tion (event) to stop integration. While integrating from 
the corner S mi „ to the belly <5 max , one may also compute 
the period of the solution, T — T(C), and so, invoke a 
secondary solver to solve for a S mln such that the com- 
mensurability condition T = n/N is satisfied. Such solu- 
tions have no extreme spikes as with n 3 = and instead 
yield rounder corners. The periodicity T of the solution 
does not any longer depend monotonically on the corner 
width of the roller, £ min , which results in a more com- 
plex solution space (some of the solutions look like the 
clover shapes reported by Bush et al. [17]). Interestingly, 
larger overlap of the parameter regions corresponding to 
polygons with different symmetry number N is seen, in 
particular for lower symmetries as low as N ~ 2. Details 
are discussed in [13, 14]. 



IV. TEMPORAL STABILITY OF THE 
CIRCULAR STATE 



We now introduce time dependence to the equations. 
First, we consider time-dependent version of the continu- 
ity equation (5). 



d_ 

dt 



R 2 



-Ahd9 



9_ 

2tt 



d9 + qg(9) 



[q r (9)d8 + q e (e + d9)}. (79) 



The left hand side is an estimated rate of change in the 
volume of the roller slice. The right hand side describes 
flux into and out of the volume as described in Sec. II B. 



This equation becomes in dimensionless form: 



T 



6 1- 



Now, the roller thickness S(t, 9) also depends on t, and 
the characteristic time scale T is 



T ; 



R 2 Ah _ R 2 h 2 Q 



(81) 



i.e., time to fill the whole disk with radius R and height 
Ah ~ h - Using this together with the radial force bal- 
ance equation (12) we get 



T 



8_ 
dt 



5 1- 



1 

2t 



Hi 

S 



89 



(82) 



Next, we introduce time dependence in the azimuthal 
momentum balance (47). This becomes 



n 4 dt 



_d5 

~ ~89 

+ 3aBo- 2 S 



n 2 u 3 (2-s) 

S 2 2 
p' x {S)dS 8 3 8 



£0 



where 



n 4 



p 2 89 d9 3 

gR 2 hl 
Q 2 ' 



(83) 



(84) 



We now linearize around the circular solution as before: 
Setting 8 — 5 + eS\{t, 9) and £g = 0+e£i(t, 6), and taking 
the first order terms in e, we find 



and 



n 4 at 



T(l - S ) 



II2 

*o 2 



dSi 
~dt 



1 



2nS, 



■Si 



oil 

89 



(85) 



n, 1 



:So 



(l - 3a£Bo 



- 2 \ dSi 
d9 



3a<5oBo 



89 3 



(86) 



instead of their stationary version (49) and (50). 
We now Fourier transform by 



'Sx(t,e)' 




z 








X 


exp 



St 



(1 - So)T 



ik9 



(87) 



where, due to the 27r-periodicity in 9, the wavenumber k 
has to be an integer. We find 

sx = ~C 2 x - i{C x k + C 3 k 3 )z 
sz = —ikx + C±z, 



where 



Ci 
C 2 

c 3 
c 4 



(1 - So)H 4 (1 - 3aBBo 

"n 2 „ L 1 



(1 - S Q )Tl 4 



n. 



1 



s 



6o(l - £ )3an 4 Bo- 2 
1 



2tiV 



(89) 

(90) 
(91) 
(92) 
(93) 
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where B is defined in (43). The characteristic equation 



s 2 + (C 2 - C 4 )s + (Cifc 2 + C 3 k 4 ) - C 2 Ci = (94) 

with solution 
1 

S= 2 



d - C 2 ± V(C4 + C 2 ) 2 - 4(Cifc 2 + C 3 fc 4 ) . 

(95) 

For = the solutions are s + (0) = C4 and s_ = — C 2 . 
Since C4 > 0, the spatially homogeneous (circular) mode 
is unstable, leading to an increase or decrease of S away 
from 5q. This is clearly unphysical (since our circular 
state is uniquely defined from the outset) and reflects 
the fact that flux conservation is not built into our time- 
dependent model, since (82) does not directly respect flux 
conservation. Thus, we must require that the constraint 
(16) still holds, whereby the total mass flux in and out 
of the roller is in balance. To linear order, this condition 
implies that 



<Ji(0)d0 = O. 



(96) 



In 



Thus, the mode k = must vanish and is excluded 
the following we only consider integer values k > 1. 

Let us take a look at the dispersion relations in (95). 
All the Il's and the Bond number are positive, so all the 
C's are guaranteed to be positive except C\. As dis- 
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FIG. 13. (Color online) There is no maximally unstable wave 
number k* when C\ > 0, see Eq. (98). The stability curves 
s+(k) (blue dotted) and s_(fc) (red dashed) from (95) for 
Ci = 2 (above) and d = -2 (below) with C 2 = 2.5, C 3 = 0.1 
and C4 = 0.3. The dashed lines show the curves for contin- 
uous k, but only integral values of k, corresponding to the 
crosses, are allowed due to 2-7r-periodicity in 8. The mode 
k — is treated separately due to the global conservation law 
(96). 

cussed in Sec. II D 4, surface shape parameter B is pos- 
itive when Ri grows with 5 - this is the case at least 
close to the circular jump where 8 ~ 5q. Then, the sign 



of C\ only depends on the relative magnitude of B and 
Bo 2 : when 3aB > Bo 2 , we have C\ < 0, and surface 
tension is "active" in the sense that it then can cause a 
Rayleigh-Plateau-like instability. This instability crite- 
rion is equivalent to 



Rh 

IT 



< W. 



(97) 



At present, we do not have much data to check this cri- 
terion, but from the case presented in FIG. 3 we can esti- 
mate the two curvature radii and thereby the magnitude 
of B « 27 (see Appendix A). Note however, that since 
the values entering B should really be averages over the 
surface, we somewhat overestimate the curvatures and 
therefore likely underestimate B. Ethylene glycol has a 
surface tension of a rj 50 x 10~ 3 N/m and a density of 
1100 kg m~ 3 ; these numbers yield a capillary length of 
l c Ri 2.2 mm and a Bond number Bo w 2.3. In effect, 
the left hand side of (97) thus becomes ~ 32, whereas 
the right hand side is 81. This demonstrates that the 
circular state should indeed be unstable towards the for- 
mation of polygons. The flow-rate is Q ~ 40 ml s~ x 
and using the parameters listed in Appendix A, we have 
n 4 = gR 2 hlQ~ 2 ri 0.69; this yields a positive prefactor 
in Ci of (1 - S )U 4 r; 0.4. Although the two sides of (97) 
appear similar, these estimates show that it is plausible 
that C\ may take on negative values. 
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FIG. 14. (Color online) The stability curve s+(fc) (dashed 
line) from Eq. (95) shown for stability coefficients Ci « 
-0.603, C* 2 w 0.415, C* 3 » 0.0156, and C* 4 Ri 0.374, based 
on experimental parameters and the height profile of a pen- 
tagon (see Appendix A). Note that the constant branch 
s(k > k ) = (C 4 - C 2 )/2 Ri -0.02 is negative. We then find 
that the maximally stable wave number (Eq. 98) is k* ~ 4.4; 
thus the most likely symmetry to occur has symmetry k = 4 or 
k — 5 which matches well the expectation. The crosses mark 
the values of s for integer k, respecting the 27r-periodicity. 
The mode k — is treated separately according to the global 
conservation law (96). 

Let us look at the dispersion relations in more de- 
tail. First we note that s(k) has reflective symmetry 
around k = 0. In Figs. 13, 14 we therefore only dis- 
play the positive fc-axis. Along this positive /c-axis, 
there are (at most) three distinct values that occur for 
the dispersion relation: (i) the maximal value of s(k) 
at k* , defining the two neighboring modes of strongest 
growth (integer values), (ii) the root of the dispersion 
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relation s(fco) = 0, and (iii) the point k c where the two 
branches s_(fc c ) = s+(fc c ) meet. For k > k c , the discrim- 
inant D = (C 4 + C 2 ) 2 - 4(Cifc 2 + C 3 k 4 ) becomes nega- 
tive and the real part assumes simply a constant value, 
Re(s) = (C4 — C2V2, which defines a further symmetry 
axis for the dispersion relation. Note that, by definition, 
k is identical to Eq. (53). 

For the unstable case {C\ < 0) we expect the symme- 
try breaking predominantly to occur at the maximally 
unstable wave number k* . Re(s+) is maximal when the 
term |Ci|fc 2 — C3A: 4 is maximal, and so 



k* = 



-Ci 
2C, 



baB- Bo 2 
GcySq 



(98) 



Clearly, such a maximal value only exists when C\ is 
negative, i.e. when the condition 3aB > Bo 2 in (97) 
is satisfied, as illustrated in Fig. 13. For small Bond 
number, or strong instability, (98) simplifies to 



(99) 




and if the roller cross section is close to cylindrical, we 
find k* ~ ^o" 1 ) i- e - a wavelength proportional to the 
width of the roller as in the Rayleigh-Plateau instability. 
This fits nicely with the observed bifurcation sequence, 
since, experimentally, polygons with many corners (up 
to 13) are seen close to the transition, where the roller is 
thin, and lower order polygons only occur later when the 
width is larger. 

In FIG. (14), we show the dispersion relation using 
parameter estimates from Appendix A which are based 
on the height profile of the pentagon shown in FIG. 3. For 
this case, we obtain k* 4.4 and fc ~ 6.23. With this 
value of k* , it is most likely that polygons with either 
k = 4 or k = 5 corners occur, in accordance with the 
experimental observation. 



V. DISCUSSION 

We hope that we have been able to convince the reader 
that the type II hydraulic jumps and their intriguing 
polygon offspring can be understood within the phe- 
nomenological framework presented here. We have given 
a detailed presentation of the model and the underlying 
assumptions about the flow. In particular, we assume 
that the inner edge of a "roller" structure inside the jump 
describes the shape of the type II jump. In the circular 
type II regime, the internal flow of this roller structure is 
similar to a toroidal vortex. The local width of this roller 
is one of three basic variables in our model. The two other 
variables in our model are the radial and tangential flow 
rates across the jump and inside the roller, respectively. 
A polygonal shape is obtained when the width loses cir- 
cular symmetry and becomes a periodic function of the 
polar angle, while the outer edge of the roller remains 



circular, as observed in experiments. Also, for the polyg- 
onal states, a net flux is carried inside the roller from the 
sides to the corners. The model is formulated by consid- 
ering a radial and a tangential force balance across the 
roller structure and the polygon shapes occur primarily 
by competition between viscous and gravitational effects. 
We have shown that one variant of the model is exactly 
solvable, even in the strongly non-linear limit, and that it 
has polygonal solutions which are quite similar to those 
observed experimentally. The resulting bifurcation di- 
agram is qualitatively similar to experimental measure- 
ments. 

The question of stability has been addressed as well. 
Here, surface tension plays an important role. For small 
deformations around the circular state, we can include 
surface tension by considering the Laplace pressure inside 
the roller. This pressure can be estimated by introduc- 
ing a parameter that characterizes how the curvature of 
the roller changes with deformation. When this parame- 
ter is positive, which is expected for the circular type II 
hydraulic jump, surface tension effects will tend to desta- 
bilize the circular state if the length scales at the jump 
(e.g. jump height h Q and radius R) are of the order or 
smaller than the capillary length. Thus a polygon state 
will emerge with a wave length of the order of the width 
of the roller, in analogy with the Rayleigh-Plateau insta- 
bility. 

To include surface tension beyond the linear approxi- 
mation would require more detailed knowledge about the 
height profile of the polygon states. Hopefully such re- 
sults will soon be available — indeed, one of the aims of 
the present paper is to stimulate renewed research in this 
direction. 

In the present study, we have neglected the rotational 
motion inside the roller which is clearly observed in ex- 
periments. The stability of a rotating liquid column with 
free cylindrical surface has been investigated in other 
studies [18, 19]. The relative importance of the rotational 
kinetic energy to the surface energy is characterized by 
the Weber number We$i = pa 3 £l 2 /a, where f2 is the char- 
acteristic rotational frequency and 2a — R 5 is the width 
of the roller. For an infinitely long cylinder, subject to ax- 
isymmetric disturbances (neglecting disturbances in e.g. 
azimuthal direction) , a necessary and sufficient condition 
for stability is given by k 2 > 1 + Wen [18], where k is 
our wave number (this result is valid for potential vortex 
flows and solid body rotation), and thus rotation could 
widen the spectrum of unstable modes. Compared to the 
roller in the hydraulic jump, these vortices are of course 
strongly idealized, and we do not know at present how 
to include them in a convincing way. 
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Appendix A: Typical Parameters and Coefficients 



is: Q w 40 x 1CT 6 m 3 /s, h, 
m, v w 10 -5 m 2 /s, g w 9 



A 'typical' set of parameters used in the experiments 

s 5 x 1CT 3 m, hi^5x 1Q- 4 
m/s 2 , a w 50 x 10~ 3 N/m, 
and p s» 1100 kg/m 3 . Clearly, the assumption /ij -C /i 
used throughout the article is then valid. With the above 
parameters, we have a Bond number of Bo = 2.3 (Bo 2 = 
5.4) and l c ps 2.2 mm. 

Furthermore, using the values from the height profile 
in Fig. 3, we read off the following numbers: the outer 
edge of the roller R ps 3 x 10~ 2 m; the width of the 
roller between two corners (belly) RSf, — 14 mm, and 
in the corner RS C = 11.5 mm; and the corresponding 
curvature radii Rb = 4.5 mm and R c — 2.5 mm. Using 
the approximations So = (Sb + S c ) /2 s» 0.426 and 



B = 5 



1 



dS pi(5) 
lSb + 5, 



So 



2S b - S c \pi(5 c ) pi(5b) 



we estimate that B ps 27. Note that the values entering 
B should really be averages over the surface; therefore 
we somewhat overestimate the curvatures and therefore 
likely underestimate B. 

The scaling factors ci,C2,C3 and C4 are unknown, but 
using the above experimental parameters and the esti- 
mates from the height profile, we can at least estimate 
some of them. Using the scaling R = ciQ 5 ' 8 ^ -3 ' 8 *? -1 ' 8 
proposed in [15] for the outer edge of the roller, R, we es- 
timate that C4 ~ 0.3; moreover, we can match the defini- 
tion Sq = 27rni with c\ ~ 1.13. The constants c%, C3 origi- 
nate from the linearization of the frictional force (26), and 
it is reasonable to assume that they are of order unity, 
too. With these estimates, the non-dimensional num- 
bers become: II 1 w 0.0678, U 2 « 0.00363, n 3 ps 1.306, 
n 4 w 0.689 , 5 « 0.426, <p « 11.6. The stability coeffi- 
cients are then C x « -0.603, C 2 ~ 0.415, C 3 « 0.0156, 
and C4 s» 0.374. With these estimated coefficients, we 
find the spectrum for stable wave numbers as shown in 
Fig. 14, with a maximally unstable wave number k* ps 4.4 
and s+ has a root kg ps 6.23. Evaluating the wave num- 
ber of the linearly perturbed circular jump in Eq. (53), 
we have k ps 6.23 which is at least quite close to a pen- 
tagon. However, we note again that these numbers can 
only serve as rough estimates due to the lack of experi- 
mental data (i.e. to determine c\, 02,03, C4 and B). 



Typical frequencies for the rotational motion in our ex- 
periment are of the order of 5 Hz, so £1 ~ 5 x 2tt rad/s. 
The roller structure has a width of order RS /2 « 6 mm 
and a height of h /2 = 2.5 mm. Using for the aver- 
age roller radius an estimate of (a — h a + R5q)/A w 4.5 
mm we arrive at an estimated Weber number of Wen = 
pa 3 il 2 /a ~ 1.9; notice though that the Weber number 
scales with the cube of the average roller width a and is 
thus very sensitive to errors in the estimate of a. It is not 
clear that the vorticity extends all the way out to where 
the surface attains its maximal value (this is how we es- 
timated RSo), and therefore we may be overestimating 
a. 

We define standard dimensionless numbers as follows. 
Defining the velocity before the jump as u- — — — 
get the Reynolds number 



2-irRhi 



WO 



Rc = 



Qh Q 
litRhiV 



212, 



(Al) 



where we take h a as the characteristic scale, similar to 
RS. The Froude numbers are 



Fri 



Q 



I ghi 2irRhi\/ghi ' 



so 



Fr? = 



37.037. 



1 {2nYgR^ 
Similarly, we have the downstream Froude number 



O 2 

Fr 2 = ^ „„. „ w 0.037. 



{2*YgR 2 hl 



(A2) 



(A3) 



(A4) 



These dimensionless numbers are evaluated using the 
above choice of parameter, with R = 30 mm for the 
height profile in Fig. 3. With these definitions, (f> from 
Eq. (76) is 



1 



R h o 1 

The Weber number is defined as 



Frr 2 Fr" 2 Re 2 . 



We = 



a/h 



(A5) 



(A6) 



because h Q is similar to RS, which is similar to radii of 
curvature. Thus 



We 



ph a Q 2 



{2n) 2 R 2 h 2 a 



(A7) 



Appendix B: Nomenclature 

The variables and dimensionless constants used 
throughout the paper are summarized in Table I. 
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Symbol 


Meaning 


a 
P 
P 

V 

Q 


surface tension 

fluid density 

dynamic viscosity 

kinematic viscosity 

volumetric flow rate from nozzle 


(r,z,9) 
u = u(r, 9, z) 
h = h(r,z,8) 
q r 

qe 

Cr = q r /Q 
£e = qe/Q 

rM 

5(9) ~ R ~ r ^ e) 


cylindrical coordinates 

flow velocity 

fluid surface height 

flow rate in radial direction 

flow rate in azimuthal direction 

non-dimensional radial flow rate 

non-dimensional azimuthal flow rate 

jump radius 

non-dimensional jump width 


hi 
h 

Ah = h — hi 

So 

R 

a = h / R 


inner fluid height 

outer fluid height 

height difference across jump 

non-dimensional (circular) jump width 

outer radius of roller structure 

aspect ratio of jump 


dA f 

dA b (9) 

dA r 

At{9) 

Cl , C2 , C3 , C4 

Rl 

R 2 

pi = Ri/R 

B = S (d/dS(l/p 1 ))s=s 


free surface area of control volume 

with infinitesimal angle d9 (inf. CV) 

bottom area of inf. CV 

back side of inf. CV 

lateral area of inf. CV 

geometrical constants 

principal radius of curvature (r,z)-plane 

principal radius of curvature (« (r,#)-plane) 

non-dimensional radius of curvature 

surface shape parameter 


dF r 
dFg 


radial force per infinitesimal angle d9 
azimuthal force per infinitesimal angle d9 


N 
k 


polygon symmetry number 
wave number 


n x 

n 2 ,n 3 
n 4 

Bo = ho/lc 

l c = yJo/(pg) 


hydrostatic vs. shear force (radial) 
shear vs. hydrostatic forces (azimuthal) 


Bond number 
capillary length 


T 

Ci, C2, Cz, Ca 


filling time: disk of radius R with height Ah 
parameters used in stability analysis 


x — 8 /So 
y = 27r^ 

4> = (2^) 2 n? /2 n 2 - 1/2 
c,v 

Q = 9/<)) 


jump width (solvable model) 
transverse flow rate 
"mass" 4> 2 

energy, potential function of " mass " <j) 2 
rescaled angle 



TABLE I. Variables and parameters used in the model. 
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